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■ Abstract. Cross-correlation function (CCF) has become the standard 

tool for extraction of radial- velocity and broadening information from 
high resolution spectra. It permits integration of information which is 
common to many spectral lines into one function which is easy to cal- 
culate, visualize and interpret. However, CCF is not the best tool for 
many applications where it should be replaced by the proper broadening 
function (BF). Typical applications requiring use of the BF's rather than 
CCF's involve finding locations of star spots, studies of projected shapes 
J> ' of highly distorted stars such as contact binaries (as no assumptions can 

be made about BF symmetry or even continuity) and [Fe/H] metallicity 
determinations (good baselines and avoidance of negative lobes are essen- 
tial). It is stressed that the CCF's are not broadening functions. The note 
' concentrates on the advantages of determining the BF's through the pro- 

cess of linear inversion, preferably accomplished using the Singular Value 
Decomposition (SVD). Some basic examples of numerical operations are 
^ ■ given in the IDL programming language. 

D.' 

6 

52 ■ 1. Convolution and cross-correlation 

03 

Convolution is an operation that the nature does for us. We seldom see "naked" 
k>( | functions. These could be a convolution of a spectrum with the spectrograph's 

j_j ■ instrumental profile or with the radial component of the micro-turbulence ve- 

locity field in the stellar atmosphere or with a broadening function due to rapid 
rotation of a star. Thus, instead of a function f(u), we observe a function h{x) 
which is a convolution with some other broadening function (BF), g{x): 

/+oo 
f(u) g(x -u) du = f{x) * gix) 
-co 

This natural process can be easily simulated in numerical packages (examples in 
the IDL programming language are marked by a command-line prompt IDL>) 
either by a special operator: 
IDL> h = convoKf , h) 

or through the Fourier-transform multiplication and its inverse: 
IDL> h = float (fft(fft(f ,-l)*fft(g,-l) , + D) 

Cross-correlation is an operation which for real functions differs from the 
convolution really only in the symmetry of the arguments. For complex functions 
things are slightly different (real and imaginary parts have different symmetries) , 
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but astronomers observe real spectra so we do not have to worry about the 
mathematical nuances. 

/+oo 
f(u) g{u + x) du = f(x) -k g{x) 
-oo 

The cross-correlation (note the different asterisk above) function can be com- 
puted numerically using: 
IDL> lag = findgen(201) - 100 
IDL> c = c_correlate (f ,g, lag) 

2. Broadening functions 

Suppose we observe a sharp-line (S(X)) and a broad-line (-P(A)) spectra and 
we want to determine the broadening and any other differences which make the 
latter spectrum more interesting than the former. The function B(X) can be 
a rotational broadening function for a single, rapidly-rotating star, or a more 
complex profile for two components of a binary, or for a star with spots (where 
they would show as indentations in the function). The sharp- line spectrum S(X) 
is not free of some broadening. This can be the thermal broadening of lines or 
micro-turbulence or some other mechanisms; we call them jointly T(A). Thus, 
schematically, the sharp-line spectrum can be written as: 

S(\) = (£a i 8(\ i ))*T(\) 

i 

while the broad- line spectrum, broadened additionally by B(X), can be written 
as: 

P(A) = S(X) * B(X) = (J2 ai5{Xi)) * T(A) * B(X) 

i 

The cross-correlation (CCF) with the sharp-line spectrum is frequently 
taken as an estimate of B(X): 

C(X) = S(X)*P(X) = S{X) * {S{X) * B{X)) = T(A) * T(A) * B{X) 

The new function £>(A), is not identical to B(X), because it inherits the common 
broadening components (such as thermal, micro-turbulence, instrumental) from 
both spectra. Tonry & Davis (1979) showed that if those common components 
are represented by Gaussians, the addition is quadratical, which for these func- 
tions means repeated convolutions. Thus, CCF cannot be really used to replace 
the broadening function. But it can give us some approximation of it and will 
remain a useful tool to have some preliminary estimate on the degree of the line 
broadening. For symmetrical broadening functions, it will remain the simplest 
tool to determine the radial velocities simultaneously from many spectral lines. 

The differences between the BF and the CCF can be seen when an artificial 
broadened spectrum is created by a convolution and then the resulting spectrum 
is subject to the CCF operation. The result (Figure 1) is obviously different from 
the BF: The CCF shows negative baseline excursions and, most worryingly, it 
shows the "peak-pulling" effect which would lead to an under-estimate of the 
individual component velocities. While this last problem can be overcome by 
applying the TODCOR technique (Zucker & Mazeh 1994), we clearly see that 
the CCF is not the BF. 
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Figure 1. The left panel shows a sharp- line spectrum (upper part) 
and the result of convolving it with the broadening function for a con- 
tact binary (lower part); this function is shown in the right panel by 
the dotted line (BF). The cross correlation of the broadened spectrum 
(with added noise to have S/N = 100) with the sharp spectrum is 
shown in the right panel by the continuous line. Notice the negative 
baseline and the peak-pulling in the CCF. 



3. The Fourier transform de-convolution 

Some attempts to determine the broadening functions (Anderson et al. 1983) 
utilized the well known property of the Fourier transforms of a correspondence 
between convolutions and multiplications in the two relevant domains. Thus, a 
convolution: 

P(\) = S(X) * B(X) 

transformed with the Fourier transform T changes into a product of the trans- 
forms: 

T {P(A)} = T {S(X)} ■ T -{B(X)} 
Therefore, the broadening function can be restored with: 

B(\)c^{T{P(\)}/F{S(X)}} 

This can be done in a compact way with: 

IDL> b = float(fft(fft(p,-l)/fft(s,-l) , + D) 

While the mathematical background is simple and easy, the practice is 
just the opposite. First, the resulting B(X) spans the whole spectral window, 
so that one determines a lot of zeroes; there is no "compression" information 
whatsoever. But, more importantly, the division operation usually ... does not 
work: the high frequency noise becomes amplified and some sort of frequency 
filtering is needed. The result may then actually depend on the applied filter. 

Some authors (see for example the large collection of works of D. F. Grey) 
use spectra transformed into the frequency domain without the division step. 
This can be done for broadening mechanisms describable by simple functions or 
obeying some symmetries, but fails in cases of spots or of BF's of close binary 
systems. 
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4. Convolution in the formalism of linear equations 

There are two main issues that re-casting convolution into a set of linear equa- 
tions can resolve. These are: (1) How to channel information over the whole 
spectrum (say 2000 pixel long) into the BF window (say 200 pixels long)? 
(2) How to utilize all information contained in sharp-line spectra and remove 
the influence of the noise in the continuum (which carries no information)? 

The convolution can be written as an over-determined system of linear 
equations which link a sharp-line spectrum S(n), via the broadening function 
B(m), with the broadened spectrum P(m). The mapping is through the "design 
matrix" Des(m,n) which is formed from the sharp line spectrum S(n) by con- 
secutive vertical shifts by one element. In IDL, this can be done with a simple 
routine: 

function map4,s,m 

; m - must be odd, n must be even 

n = n_elements (s) & t = fltarr(m) # f ltarr (n-m+1) 

; t(m,n-m) = t (small , large-small) dimensions 

for j = 0,m-l do for i = m/2,n-m/2-l do t ( j , i-m/2)=s (i-j+m/2) 

return, t 

end 

An example of using this routine to create a design matrix for a 201-pixel 
long window would be: IDL> des = map4(s,201). The program spectra must 
accordingly be trimmed to n-m+1 with: IDL> p = p(m/2:n-m/2-l). The sys- 
tem of equations has a familiar form of the over-determined linear set: 
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5. Singular value decomposition (SVD) 

One of the traditional ways of solving the system of equations above would be 
to transform it into a system of the "normal" equations of the size reduced from 
m x (n — m) to m x m by multiplication of both sides by the transpose of the 
design matrix. The result would be the BF defined in the least-squares sense, and 
it is possible to stop at this point. However, the Singular Value Decomposition 
technique also gives us such an answer, but - in addition - makes it possible to 
remove the influence of the continuum and its noise. 

The SVD technique is beautifully described in the "numerical techniques 
Bible" of Press et al. (1986). They present it as a somewhat magic black box 
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and for most users it is just fine. If you want to learn how the technique really 
works, then the books of Golub & Van Loan (1989) or Craig & Brown (1986) 
are probably the best references. 

The essence of the SVD is the property that one can represent any matrix 
by a product of 3 matrices; in our case: Des = UWV T . These matrices are: 
the column ortho-normal U and V and the diagonal matrix W (this is really a 
vector containing the diagonal elements). The property of the columns in U and 
V T is that the following products, U T U = I and V T V = I, give the unity array 
1 (1 on the diagonal). 
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In IDL, the operation is represented by: IDL> svdc, des, w,u,v, /double. Here, 
des is the only input quantity and the remaining parameters are what the routine 
produces as output. The keyword /double is for higher precision and is optional. 

One can check the correctness of the operations by the following commands: 
IDL> wf = fltarr(m,m) 
IDL> for i = 0,m-l do wf(i,i) = w(i) 
IDL> des_check = u ## wf ## transpose (v) 

The three new matrices are all invertible: U and V are ortho-normal arrays, 
so that their inverses are just transposes, while the diagonal array W is replaced 
by a diagonal array Wl, with the diagonal elements containing the inverses, 
Wlij = l/wf. 
IDL> wl = fltarr(m,m) 

IDL> for i = 0,m-l do wl(i,i) = l./w(i) 

The solution is given by: B = VW\{U T P) or schematically: 
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Figure 2. The left panel shows the run of the singular values for the 
sharp- line spectrum in Figure 1. The right panel gives the solutions 
for the broad-line spectrum in the same figure (with added noise at 
S/N = 100) utilizing the first 5, 10, 20 and 30 singular values. 



Numerically: IDL> b = ref orm(v##wl## (transpose (u)##p) ) . A routine 
makes this simpler: IDL> b = svsol(u,w,v,p, /double). The elements of B 
are all independent, so that any - even strange or discontinuous - broadening 
functions can be restored as no condition imposed on the smoothness or sym- 
metry of the result. Note that, if only one sharp-line template is used, the 
decomposition operation svdc is done only once, for possibly many broad-line 
spectra p, each giving a separate solution b. 



6. Advantages and disadvantages of the SVD approach 

On the positive side: (1) The problem can be treated as a linear-equations. 
(2) An "inverse" of the rectangular array Des is possible. (3) The solution of 
B is defined in the least-squares sense (shortest modulus). (4) The result is the 
real broadening function. But there are also minuses: (5) One must solve a large 
system of, say, 2000 equations for 200 unknowns. (6) One must know a priori 
how many unknowns. (7) Initially, the results may turn out quite poor, because 
of the presence of plenty of linearly-dependent equations in the system (parts of 
spectra where the featureless continuum provides no broadening information). 

The SVD approach offers a simple resolution of (7) as it permits removal 
of the effects of the continuum in an objective way. The key element here are 
the singular values contained in the diagonal of W. Since the solution involves 
l/wi, small values in Wi spoil the solution. These are exactly those problematic 
values that one wants to avoid. Thus, by rejecting of small values of Wi, one can 
(i) remove the linearly dependent equations, (ii) diminish the influence of the 
noise from the continuum, (iii) reduce the influence of the computer round-off 
errors (which enter multiplied by the order of the problem) and (iv) reduce the 
number of the unknowns (because the system is usually not over-determined at 
all). All these properties are related to the "conditioning" of the array Des. 
The reader is directed to the source texts on this subject for further reading. 

The important factor is m&x(wi) / mm(wi) which provides an estimate on 
how many of the singular should be used. In practice, one can keep on adding 
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more wi and see the successively better solutions. The diagonal arrays W 
and WI will have then elements: Wi = wq, wi, W2, w^, w m -\ and wli = 
1/wo, 1/wi, 1/wk-i, 0, with k (we call it the order of solution) spanning 
the whole range to m — 1. In IDL, this can be done by forming a square 
matrix of solutions: 
b = f ltarr (m,m) 
for i = 0,m-l do begin 
wb = f ltarr (m) 

wb(0:i) = w(0:i) ; first i+1 singular values used, rest zero 
b(*,i) = svsol(u,wb,v,p, /double) 
end for 

We note in passing, that the arrays U and V have very special properties 
as they contain the basis vectors in the spaces of the spectra and broadening 
functions, respectively. One can see this by analyzing a diagonalized system: 
W Z = D, obtained by keeping the same W, with Z = V T B and D = U T P. 
The solution of the diagonalized system would be then: Z = D/W, but, in 
practice, the diagonal of W is a vector, so that Zi = dj/iUj. Plotting the columns 
of U and V can tell one a lot about the conditioning of the solution. 

7. A few notes on the SVD solutions 

The first question is: How far in k should one go and where to stop? The essential 
operation is to plot (usually in log units) the vector W (Figure 2). There are 
3 parts of it: (1) the good, large singular values, (2) the small values usually 
representing the noise in the spectrum 5(A) and (3) the numerical errors. You 
may want to stop no later than at the kink below the good part. But the real 
"quality control" is the fit to P. If the error of the fit stops decreasing, you have 
found the right point (Figure 3). Beyond that point, you will start fitting the 
noise! However, it is useful to analyze the solutions for different orders k and see 
how they first improve and then get worse. Sometimes the fit will remain poor, 
in spite of the leveling of the error curve; this usually means a wrong choice of 
the sharp-line spectrum. The standard error of the fit can be calculated from: 
sig = f ltarr (m) ; error 

pred = des ## transpose (b) ; predicted fits 

for i=0,m-l do sig(i) = sqrt (total( (pred(i , *) -p) ~2) /m) 

Usually, even very low order solutions are well defined (Figure 2), but stop- 
ping early is not always advisable because this leads to a loss of resolution, as 
then not all basis vectors contribute. The solution which - in principle - has all 
elements in B independent suddenly acquires inter-element correlations. Thus, 
it may be advantageous to go to a highest possible order (A; = m) and thus insure 
that elements of B are uncorrelated, and then decrease the noise by smoothing. 

One should be aware that the errors may be under-estimated for the case 
of truncated (k < m) solutions. While the prescriptions of Rix & White (1992) 
and Rucinski, Lu & Shi (1993) are based on the theory of the full SVD, the 
error analysis for the truncated case has not yet been done. In this situation, it 
may be advantageous to utilize techniques of the external estimates, such as the 
bootstrap or Monte Carlo. This subject certainly requires more work ... 
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Figure 3. The left panel shows the mean standard error of the fit 
for our example solutions utilizing progressively more singular values. 
It is obvious that the first 10 - 15 singular values give an adequate fit 
(compare with Figure 2). The right panel shows the solution utilizing 
all 200 singular values, convolved with a Gaussian. This approach may 
offer a better control over the final resolution of the BF. 

8. Conclusions 

One can position the SVD technique of linear broadening function restoration 
as located between the cross-correlation and the direct modeling of the spectra. 
The BF's determined with it are much better defined than the CCF's, and are 
true broadening functions, not their proxies. They also integrate the geometrical 
information from a spectrum, but give well defined baselines without the CCF's 
negative fringes. They do require a bit more computer work, but several numer- 
ical packages can easily handle large linear systems of equations involved here. 
Obviously, they cannot replace the spectrum synthesis, if these are needed, but 
can be a useful tool in their preparation. 
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